package optimization;

import static optimization.LinearAlgebra.*;

import java.util.ArrayList;
import java.util.List;


import common.statistics.Series;
import java.util.logging.Level;
import java.util.logging.Logger;

public class RetrospectiveApproximation {

	StochasticFunction f;
	VectorFunction cons;


	List<Integer> nSamples = new ArrayList<Integer>();
	List<double[]> approxRoots = new ArrayList<double[]>();
	
	public int iter=0, N;  // iteration, current # of samples
	public long randSeed = 10000000L;

	public double hscale, term;  // h scaling,  termination threshold
	public double[] xa, h; // approximate root for iteration, delta vector


	
	public RetrospectiveApproximation(StochasticFunction f, VectorFunction cons, double[] x0, double[] h) 
	{
		this.f = f; 
		this.cons = cons;
		this.xa = x0.clone();
		this.h = h;
	}
	
	public double[] findApproxZero()
	{
		NelderMeadS9RS nm = new NelderMeadS9RS(f, cons, term);
		nm.setStochasticParameters(N, randSeed);
		return nm.search(xa, vMul(hscale, h));
	}

	/**
	 * Override to modify hscale updating
	 * @param iter
	 * @return 2 / 2^(i*0.7)
	 */
	protected double hscale(int iter)
	{
		return 2 / Math.pow(2, iter*0.7);
	}
	
	/**
	 * Override to modify termination threshold
	 * @param iter
	 * @return hscale(iter) / 20
	 */
	protected double term(int iter)
	{
		return hscale(iter)/20;
	}

	/**
	 * Override to modify sample size generation function
	 * @param iter
	 * @return 16 + 2^i
	 */
	protected int sampleSize(int iter)
	{
		return (int)(16*Math.pow(2,iter)); 
	}

	
	public void beginIteration()
	{
		// set new random seed
		randSeed++;
		
		// determine sample size
		N = sampleSize(iter);

		// determine h, term for iteration
		hscale = hscale(iter);
		term = term(iter);		
	}

	public void endIteration()
	{
		// save approx root
		nSamples.add(N);
		approxRoots.add(xa);
		
		iter++;
		// use approx root as next initial guess
//		x = xa;
	}
	
	/**
	 * Search result for <iterations> iterations.
	 * @param iterations
	 * @return root estimate
	 */
	public double[] search(int iterations)
	{		
		for (int j = 0; j < iterations; j++) {
			beginIteration();
			
			// find approximate root
			xa = findApproxZero();
			
			// show stats
			try {
				Logger.getLogger("com.ristolaakso").logp(Level.INFO, Thread.currentThread().getStackTrace()[0].getClassName(), Thread.currentThread().getStackTrace()[0].getMethodName(),
						  String.format("RA%d   X %.5f, %.5f   (N=%d, rnd=%d, h=%.5f, nu=%.5f)\n",(j+1), xa[0], (xa.length > 1? xa[1]: -1.0), N, randSeed, hscale, term));
			} catch (Exception e) {
				System.out.println("Print Error: " + e.getMessage());
			}
			endIteration();
		}		
		return estimateRoot();
	}

	public double[] estimateRoot()
	{
		// IRA; result is weighted sum of approximate roots
		double[] rootEst = new double[approxRoots.get(0).length];
		double nSum = 0;
		for (int i = 0; i < nSamples.size(); i++) {
			int nS = nSamples.get(i);
			double[] approx = approxRoots.get(i);
			rootEst = vAdd(rootEst, vMul(nS, approx));
			nSum += nS;
		}
		rootEst = vMul(1/nSum, rootEst);

		Logger.getLogger("com.ristolaakso").logp(Level.INFO, Thread.currentThread().getStackTrace()[0].getClassName(), Thread.currentThread().getStackTrace()[0].getMethodName(),
			String.format("RA**   X %s\n", LinearAlgebra.vPrint(rootEst)));

		return rootEst;		
	}

	
	
	public static RetrospectiveApproximation nmEmuRA()
	{
		final double[] ystar = new double[] { 5.35, 4.41 };
		double[] init = new double[] { 1.8076, 0.0516 };

		StochasticFunction f = new StochasticFunction() {
			@Override
			public double evaluate(double[] p, int ns, long rngseed) {
				Series s = NMEmu.emulate(p[0], p[1], ns, rngseed);
				double dm = s.mean() - ystar[0];
				double ds = s.stdev() - ystar[1];
				return Math.sqrt(dm*dm+ds*ds);
			}
		};
		VectorFunction cons = new VectorFunction() {
			@Override
			public int yDim() {
				return 2;
			}
			@Override
			public double[] evaluate(double[] p) {
				if (p[0] < 1) p[0] = 1;
				if (p[0] > 6) p[0] = 6;
				if (p[1] < 0) p[1] = 0;
				if (p[1] > 1) p[1] = 1;
				return p;
			}
		};
		double[] h = new double[]{ 0.1, 0.05};
		
		RetrospectiveApproximation ra = new RetrospectiveApproximation(f, cons, init, h);
		return ra;
	}

	
	
	public static void main(String[] args) {

		RetrospectiveApproximation ra = RetrospectiveApproximation.nmEmuRA();
		double[] rootEst = ra.search(10);

		Logger.getLogger("com.ristolaakso").logp(Level.INFO, Thread.currentThread().getStackTrace()[0].getClassName(), Thread.currentThread().getStackTrace()[0].getMethodName(),
			"** FINAL ESTIMATE **");
		Series s = NMEmu.emulate(rootEst[0], rootEst[1], 1024, 0);
		Logger.getLogger("com.ristolaakso").logp(Level.INFO, Thread.currentThread().getStackTrace()[0].getClassName(), Thread.currentThread().getStackTrace()[0].getMethodName(),
			String.format("RA**   X %.5f, %.5f   TS %.5f, %.5f  (N=%d, rnd=%d)\n", rootEst[0], rootEst[1], s.mean(), s.stdev(), 1024, 0));
	
	}

}
